function [ y,x ] = RK_ode44( dy_func,y0,x0,h,x_range )
%RK_ODE44 Summary of this function goes here
%   Detailed explanation goes here
x_n=x0;
y_n=y0;
h_half=h/2;
x=[];
y=[];
while x_n<x_range(2)
    if x_n>=x_range(1)
        x=[x x_n];
        y=[y y_n];
    end % if x_n>x_range(1)
    k1=dy_func(x_n,y_n);
    k2=dy_func(x_n+h_half,y_n+h_half*k1);
    k3=dy_func(x_n+h_half,y_n+h_half*k2);
    k4=dy_func(x_n+h,y_n+h*k3);
    y_n=y_n+h/6*(k1+k2*2+k3*2+k4);
    x_n=x_n+h;
end % while x_n<x_range(2)
end % function [ y,x ] = RK_ode44( dy_func,y0,x0,h,x_range )

